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Abstract 

We present mode-coupling equations for the description of the slow dynamics 
observed in supercooled molecular liquids close to the glass transition. The 
mode-coupling theory (MCT) originally formulated to study the slow relax- 
ation in simple atomic liquids, and then extended to the analysis of liquids 
composed by linear molecules, is here generalized to systems of arbitrarily 
shaped, rigid molecules. We compare the predictions of the theory for the 
(/-vector dependence of the molecular nonergodicity parameters, calculated 
by solving numerically the molecular MCT equations in two different approx- 
imation schemes, with "exact" results calculated from a molecular dynamics 
simulation of supercooled water. The agreement between theory and simula- 
tion data supports the view that MCT succeeds in describing the dynamics 
of supercooled molecular liquids, even for network forming ones. 

PACS numbers: 61.25.Em, 64.70.Pf, 61.43.Fs, 61.20.Ja 
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I. INTRODUCTION 



In the last years, the problem of a detailed theoretical description of the dynamics of 
molecular supercooled liquids has been at the center of renewed interest. The success of 
MCT for simple liquids in describing the weakly supercooled regime has stimulated a 
considerable amount of work to extend this approach to molecular liquids. Recent contribu- 
tions include the extension of MCT to describe the rotational dynamics of one linear probe 
molecule in an atomic liquid J7|, the extension of MCT to describe the dynamics — both 
self and collective — of liquids of linear molecules [HH and the MCT extension to treat 
the dynamics of a full molecular systems using a site-site representation [ ID| . In the first 



part of this article we report one further step forward in the description of the dynamics of 
supercooled molecular liquids, i.e. the extension of MCT to describe the orientational and 
translational dynamics for liquids composed by rigid molecules of arbitrary shape. Within 
the framework of " nonlinear fluctuating hydrodynamics" this has already been achieved by 



Kawasaki ||11|| . However, the connection with the molecular correlation functions (see be- 
low), which are important also from an experimental point of view has not been worked out. 
This extension will be called molecular mode coupling theory (MMCT). 

The theoretical predictions formulated in all these approaches, and the proposed approx- 
imations are currently under investigations by several research groups. Detailed tests of all 
these approaches are requested to estimate the quality of the approximations, the predictive 
power of these different approaches (i.e. the detail of the theoretical descriptions) as well 
as the limit imposed by the complexity of the MCT equations and the feasibility of their 
solution. Furthermore, the quality of the approximations may depend on the molecular liq- 
uid under investigations, for example on the fragility of the liquid in Angell's classification 
scheme. 

In the second part of this article we present the g-dependence of molecular nonergodic- 
ity parameters calculated from the MMCT equations for a model system of a supercooled 
molecular glass-forming liquid. We then compare the theoretical predictions with equivalent 
quantities calculated from extensive molecular dynamics (MD) simulations. We choose to 
compare theory and MD results for liquid water. The choice of water as molecular liquid 
is particularly relevant for testing the quality of the MMCT approach, due to the presence 
of an extended network of hydrogen bonds and to the peculiar local tetrahedral ordering. 
The inter-molecular water- water interaction is defined by the SPC/E potential ||12|| , which 
describes the molecule as a rigid planar body and models the pair interactions as a sum of 
electrostatic and Lennard Jones terms. Due to the partial charges of the atoms the molecule 
possesses a dipole moment which is directed along the 2-fold rotational symmetry axis. 
When referring to a body-fixed frame of reference for the molecule we will always choose 
the direction of the dipole as z-axis and the x-axis in the plane which is spanned by the 
molecule. 

The test of MCT and MMCT can be done in two different ways. First, one can investigate 
the validity of some qualitative predictions which are independent of the system, i.e. which 
are the same for both theories. These are e.g. the power law dependence of the various 
relaxation times on T — T c , where T c is the critical temperature for an ideal glass transition, 
and the universal relationship between the corresponding exponents and their connection 
with the so called exponent parameter A, etc. Secondly, one may calculate the numerical 
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value of e.g. A or of the g-dependence of the critical nonergodicity parameters, critical 
amplitudes, etc. Since these quantities depend on the system one has to use MMCT in case 
of molecular liquids. In this first sense it has been shown in previous papers [P^T7| that 
mode-coupling theory appears to be the correct theoretical framework for the description of 
the a- relaxation behavior of SPC/E water. Indeed, the results from a MD-simulation can 
be satisfactory compared to the predictions of MCT. The comparison has been performed 
both for center of mass (COM) self and collective correlators JT3HIH] and g-independent 



rotational correlators ||16|| . Moreover the general predictions of MCT for the g-dependent 
molecular correlators Si n> i> n >(q, m, t) have been tested up to I — 2 in Ref. [p7| , |IS| , finding again 
a remarkable agreement, both from the qualitative and the semi-quantitative point of view. 
This wide spectra analysis has confirmed the validity of the universality of MCT predictions 
and has stimulated a deeper, more quantitative understanding of the MCT approach, which 
we present in this article. 

The structure of the article is the following: In Sec. [TJ we present the complete set of 
MMCT equations for the slow dynamics in supercooled homogeneous liquids composed by 
molecules of arbitrary shape. In Sec. [TTT] we discuss the two different approximations which 
we employ to numerically solve the MMCT equations for the general molecule case and 
discuss some of the numerical techniques used to solve the MMCT equations. Finally, in 
Sec. [TV] we solve the equations for the molecular nonergodicity parameters F [n j> n /(q,m) for 
SPC/E water and we compare the calculated predictions with the corresponding quantities 
evaluated from the analysis of the MD trajectories. 

II. THEORY 

A. Collective correlation functions 

We consider a system of N identical, rigid molecules of arbitrary shape described 
by the center of mass positions xj(t) and the orientations given by the Euler angles 
Qj(t) = (<f)j(t),9j(t),Xj(t))] j = 1,2,...,N. The microscopic density p(x, ft, t) = J2j^(^~ 
Xj(t)) 5(Q, fij(t)) in the 6iV-dimensional configuration space is expanded with respect to the 
complete set of functions given by the plane waves and the Wigner matrices D l mn (fl) to 
get the tensorial density modes 

N 

p K (q,t) = Pimn(q,t) = i\2l + l)^e iqx " (t) D l * n (^(t)) (1) 

i=i 

which are the starting point of a generalized theory. Here I runs over all positive integers 
including zero, and m as well as n takes integer values between —I and /. To simplify the 
notation we often combine l,m and n in the superindex k = (l,m,n). Our aim is to give a 
close set of equations for the matrix S: 

s.A^t) = ^(p*Mt)pAq)) (2) 

of molecular correlation functions. The strategy to derive such equations of motion is the 
same as was used for atomic liquids and linear molecules and therefore explicit calculations 
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shall not be given here. A survey of the MCT formalism for simple liquids can be found in 
1 1 and detailed calculations for a molecular system can be found in || . We restrict ourselves 
to pointing out where modifications have to be made in order to treat molecules of arbitrary 
shape. 



Already in the choice of "slow variables" for the projection operator formalism p0|j21 
one is forced to make some further considerations. In the MCT for atomic one-component 
liquids one uses the density modes p(q) and the longitudinal current density modes j(q). For 
molecular liquids this has to be modified. Besides translational degrees of freedom (TDOF) 
the molecules possess also orientational degrees of freedom (ODOF) and thus we will have 
to consider two different kinds of current densities. The translational current density modes 

N 

f^(q,t) = i\2l + 1)* E^(0^ (t) ^ n (%(i)) (3) 

i=i 

describe the change of the density p K (q,t) due to translational motion of molecules, while 
the rotational current modes 

N 

j^(q,t) = i\2l + 1)1 E<W^ W ^n + ,(^W) (4) 

are responsible for the change of the density due to molecular reorientation. Here u'j denotes 
the angular velocity in the body fixed frame. Consequently also is a vector in the body 
fixed frame. We have skipped the prime because of notational reasons. As in the theory of 
angular momentum in quantum mechanics it is more suitable to use spherical components 
/i = 0, ±1, defined by a = a z , a^ 1 = -^(^x ± ia y ), instead of the usual cartesian vector 
components a x ,a y ,a z . The connection between density and current modes is given by the 
continuity equation 

p«m= E E v?r(q)j7m (5) 

a=T,RtM=-l 



with the m-independent coefficients 



-^(q x ±iq y ) a = T,p = ±l 

q z a = T, p = 



(6) 



+ 1) - n(n + /i) a = R,p±l 
n a = R, n = 



Apart from the separation into translational and rotational currents that was already 
necessary in the case of the liquid of linear molecules we will consider all compo- 
nents of the currents as slow variables instead of the longitudinal translational cur- 
rents j^{q,t) = l/qY^u,qJn*{<t) jl^iQ^t) an d combined rotational currents j^(q,t) = 

l/^l(l + l)E,qT(q)mq,t) only. 

This increase of the number of relevant variables for the projection operator formalism 
is enforced by the following considerations. For general molecules possessing three distinct 
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moments of inertia I\ 7^ I2 7^ I3 it is an important to consider the single components 
as slow variables instead of j£ from above, because only in the first case the long time dy- 
namics becomes independent of the moments of inertia, which is demanded by experimental 
observations. It is also sensible to consider the components jJ M of the translational currents 
separately since the evolution of the density correlators is influenced by dynamic correlations 
of all components of the currents as can be seen from the following equation: 

^(p*Mt)pAo)) = -EEC(^(9)0r(?.*)jf(«)) (7) 

«a' n/j,' 

in which terms {j^*(q, tyjtf 1 (o)) occur. This dynamic coupling to the transversal 
currents is induced by the anisotropy of the molecule and exists also in the case of linear 
molecules. 

The projection operator formalism then leads to the following continued fraction for the 
Laplace transform S(<f, z) — i S(q, t)e lzt , Imz > of the molecular correlation functions: 

S(q,z) = - \zS-\q) + S~ 1 (q}K(q,z)S~ 1 (q}r 1 (8) 



K K ,M z) = £ E C®<C " (?. *) ( 9 ) 

aa' /ifi' 



k(q,z) = - [zJ-HD+F^M&zU- 1 ®] (10) 

where J^, a M (q,z) = (q)) is the matrix of static current correlations. The 

reader should note that the under-bar of k, J and M stands for the additional superscripts 
afi and a'fj,'. 

Thus the derivation of an equation of motion for S(q,t) has been reduced to finding an 
expression for the memory kernel Mfg. t) , also called fluctuating force kernel, since it is the 
correlation function of a fluctuating force. The approximation scheme of MCT consists of a 
separation of J _1 (g)M(g, t) J~ 1 (?) into a regular part m re9 (q,t) which accounts for the fast 
motion and a contribution m(<f, t) due to slow pairs of density modes. For an evaluation of 
m(q, t) the fluctuating force is projected onto pairs of density modes P« 1 (gi)p re2 (g2)- This 
projection which introduces the first approximation leads to a time-dependent four-point 
correlator and a vertex function which involves static two-, three- and four-point correla- 
tors. In a second approximation both four-point correlators, static and dynamic ones, are 
factorized into a product of two-point correlators. In a final step one can approximate the 
static three-point correlator by use of the generalized convolution approximation which is 
easily generalized from linear || to arbitrary molecules. As already pointed out above we 
will not give the complete derivation here since it is rather involved and analogous to the 
case of linear molecules. The result for m(g, t) can be expressed as follows: 

9*1:92 K2A2 

{qi,t)S K2jK > 2 {q 2 ,t). (11) 
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Here the prime denotes summation such that q\ + = Q- Besides the time-dependent 
molecular correlators there appear the number density po — N/V and the vertex functions 
which are products of 

= Y.<^M^&)^Aqi) + (i «-> 2) (12) 

where c(g*) denotes the matrix of direct correlation functions which is related to the static 
structure factors by 

<q) = — (1 - S-\q)). (13) 
Po 

The coefficients u are given by 

uZ 1K M & fi) = ( (2/l+ 2 /+ 2 j 2 + 1) ) 2 C(W 2 i; ™ im2 m) x 

Y C{lihl]n x n 2 n) a = T 

X «iniW|c(i 1 j 2 j ;ni+/injn + M ) a = jR ^ 

and Ci^^^mxm-im) denotes the usual Clebsch-Gordan coefficients [ |H| . 

The equations (j8[-|10|) together with equation (0) form a closed set of equations for the 
correlation matrix S(q,t). The vertex functions given by equations ( ]T2"| - |l4|) are determined 
by the number density p an d the static structure factors S(q), only. 



B. Tagged particle correlation function 

We will now examine the motion of a single molecule that is immersed in a molecular 
liquid. Again we have a liquid of iV — 1 identical, rigid molecules of mass M and moments of 
inertia I1J2J3 whose center of mass coordinates are denoted as Xj(t) and whose orientations 
are given by the Euler angles Qj(t). In addition we have, as the iV-th particle, a molecule 
of mass M' and moments of inertia J^J^/g. As a special case we get the self correlator 
of a selected particle within a homogeneous liquid if the tagged particle is identical to the 
molecules of the surrounding liquid. 

Besides the quantities we have already introduced in the previous subsection we have to 
consider the density modes for the tagged particle 

pl;\?,t)=e*«WD^(n N (t)) (15) 
and the corresponding tagged particle correlation function 

s&m = <p£H&t)p*m (i6) 

The slow variables for the projection operator formalism are given by the density modes 
(|T5|) and the current densities 

ji s ^(q, t) = x%(t)e^Dt n (n N (t)) (17) 
J^(q,t) = ^e^D^iSlNit)) (18) 
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of the tagged particle. Also in this case it is necessary to use all components of the rotational 
currents separately to avoid inertia effects in the long-time behavior. The results of the 
projection operator formalism are analogous to the equations for the coherent correlations. 
As a further simplification the static self-correlations are given by S^' K ,(q) = 5 KjK >. Thus we 
get 



;1 + K^(g, 



(19) 



K ( S(q,z) = EEC(^S'^>;?' aV (9>) (20) 

aa' /ifi' 



k« (q, z) = - [zJ {srl (q) + J (sr 1 (<f)M (s) (<?, z) J W ~ 1 (q)] ' (21) 

with the same coefficients q^{q) as above. 

The correlations of the tagged particle will be controlled by the coherent correlations. 
Therefore in the mode coupling approximation for the memory function m> s '(q, t), the slow 
part of {q)~M^ (<?, t) J} s ' (q) , the fluctuating forces are projected onto bilinear products 
of tagged particle and coherent density modes. In the thermodynamic limit the effect of 
the tagged particle on the surrounding liquid can be neglected and the coherent correlator 
is identical to the correlation function for the homogeneous liquid. The mode coupling 
approximation for the 4-point correlation functions thus leads to the following expression 
for the memory function 



m (s)«C."'c' 
m K, K ' 



m = (^) 2 ^E' E E <?i, g 2 )^lf (g, gi, ? 2 ) 



91,92 KljKj K2,K' 2 



>< S^qi,t)S K2A (q 2 ,t). (22) 
with the vertex functions 

v^Mqx,^) = E«SU(«T» ^ ( 23 ) 



The coefficients u are the same as given above and the direct correlation function that 
describes the interaction between the tagged particle and the surrounding liquid is defined 
by 

Po(c {s) (q)S(q))^> = (pi s) *(q)pAq)} - (24) 

In the special case that the tagged particle has the same properties as the molecules of the 
liquid (q) is just the ordinary direct correlation function of the homogeneous liquid. 

Detailed investigations for the tagged particle correlators have been done for a dumbbell 
molecule in a simple isotropic liquid 0. The equations given here are the generalization 
of this theory (for a linear molecule in a simple liquid) to the general case of an arbitrary 
shaped molecule in a molecular liquid. 
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III. APPROXIMATIONS 



The aim of our numerical investigations was to examine the long-time behavior of the 
solutions of the equations of motion presented in section |II A| , i.e. to calculate the critical 
nonergodicity parameters F(<f) = lim^ 00 S(g, t) and the transition temperature T c . As 
input for these calculations we have determined the static structure factors S(q) from a 
MD-simulation, as described in detail in a previous publication |p7 |. As discussed there, the 
C^-symmetry of the water molecule leads to the consequence that the distinct part S^ d \q) 
of the structure factors vanishes for n and/or n' odd, i.e. Si mn y m i n i(q) with n and n' odd only 
contains information about the self-correlation of the molecules. These symmetry relations 
allow for a simplification of the equations. As shown in appendix [A] the matrix equation 
splits into two parts. The matrix elements Si mnt i' m ' n '(q, t) with n and n' even are determined 
by a closed set of equations which is independent of the correlators with n and/or n! odd. 
The remaining part of the equations, i.e. for n and n' odd, is identical to the corresponding 
tagged particle equations for the self-correlators as presented in section |11B| . 

For these numerical studies it is further useful to transform the equations to the q- 
frame representation, i.e. to choose the z-axis of the laboratory frame in direction of the 
vector q. The resulting set of equations and some details of their derivation are given in 
appendix [B|. The g-frame offers the advantages that the correlation matrices depend only 
on the modulus q = \q\ and in addition are diagonal with respect to the indices m and m'. 
Thus we have to solve self-consistently a set of equations for the nonergodicity parameters 
Fin t i' n '(q,m) = Fi mn ^ mn i(qe z ). 

The main computational problem in solving the equations for the glass form factor is 
the calculation of the memory matrix rn^f M (q), due to the enormous number of terms 



in the summation of eq . ([TT|) . Of course, any attempt of numerical calculation requires the 
introduction of an upper cut-off l co in I, in order to have a finite number of coupled 
equations. The stable solution F(q, m) is found as the fixed point of the iterative solution 
of the given equations. 

It has been estimated in several MCT calculations that, in order to have a reasonable 
convergence towards the fixed point, it is necessary to perform several hundreds of iterations. 
Although it is possible to considerably reduce the number of elements in the sum of eq. ( |TTD 
by taking into account the symmetries of the molecule, a full solution for l co = 2 is still not 
feasible. 

As discussed in section || it is necessary to take into account the components of the 
rotational currents j Rfl separately to avoid inertia effects in the long-time behavior. Thus 
we have taken into account all corresponding memory functions. With respect to the trans- 
lational currents we have decided to take into account only the longitudinal components j T0 , 
i.e. all memory functions with a = T, = ±1 or a' = T, // = ±1 are neglected. 

The structure of the MMCT equations further offers the possibility of several approxima- 
tion schemes, differing in the choice of the molecular static structure factors which are taken 
into account. In this article we present MMCT calculations for two different approximations, 
which bring the numerical calculations to the frontier of the present computer facilities. In 
both approximations we neglect the third angular index n in the static quantities, i.e. we 
include as input of the calculation only the static structure factors Sioj>o(q, m) and, thus, the 
direct correlation functions Qo/o(qs to). In addition, we put Fi n ^ n i{q, m) to zero for n and n' 
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different from zero. This approximation reflects in a reduction of the number of independent 
memory kernels to be calculated, although it does not imply the simplification n = n' = 
in the memory functions. 

Examining the intensity and the temperature dependence of the static correlation func- 
tions we have given some justification for this approximation in a previous publication [|17| . 



Still it has to be noted that it is mainly motivated by the need to reduce the computational 
burden and we plan to put a significant effort in the direction of a full solution of the MMCT 
equations, including also the angular index n. 

In Ref. J17[| , it has been shown that the distinct part of the structure factors vanishes 
for odd n and/or n! . Thus, with / < 2, the approximation n = n! = essentially means 
neglecting the coupling with the correlators with |n| = 2 and/or \n'\ = 2. This approximation 
is equivalent to neglecting the third Euler angle x-> the rotations of the water molecules 
around the dipolar axes. This means that the motion of the water molecules is reduced to 
the motion of their dipole moments. Thus, we will refer to this first simplification as the 
dipole approximation. 

An even stronger approximation is defined by assuming in addition that both the static 
structure factors, the critical nonergodicity parameters and the memory functions are diag- 
onal in I and I', i.e. 

Sio,i'o(q,m) = SiQ,io{q,m)8i V (25) 
Fio,i'o(q, m) = F l0 j (q, m)5 w (26) 

m lmn,Vm'n'W} = m lmn,lm' 'n'(l)^W (27) 

Thus S(q, m) becomes a diagonal matrix with 6 non vanishing elements (the diagonal ones) 
while m(g) is still non diagonal with respect to afi and a'n'. This approximation and also 
the even stronger restriction of additional diagonality of the memory kernel with respect 
to afi has also been used in the study of a dumbbell in a simple liquid |7j and a liquid of 



diatomic molecules 22 



We have iteratively solved the equations for the nonergodicity parameters on a grid of 
100 q- vectors ranging up to llOnm -1 . Within the diagonal approximation, one iteration step 
requires about 15 minutes on one a-station. The dipole-approximation requires about two 
hours of CPU time and we estimated that a full solution of the MMCT equations including 
n would require about 4 days per iteration. 

IV. RESULTS 

We have found T^ MCT = 206 K in the diagonal-dipole approximation. At this tem- 
perature the solution of the nonergodicity parameter equations has been iterated until the 
average difference over the whole q range between the (n + l)-th iteration F^ n+1 ^(g, m) and 
the n-th F^(g, m) was of order 10~ 10 . Since the MCT-approach to e.g. a hard sphere sys- 
tem |23[], binary liquids |24| and diatomic molecules [25j has shown that MCT overestimates 



the freezing into a glassy state, we consider fortuitous the agreement between the estimated 
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rpMMCT an( j numer i ca i estimate of the critical temperature from the analysis of the 
molecular dynamics data, T C MD = 200 ± 3 K, for SPC/E water. We want to highlight that 
despite the diagonalization approximation, differently from what has been found for the liq- 
uid of Lennard- Jones dumb-bells described in Ref. [p2| , in the case of water the theory does 



not yield an unphysical splitting of rotational and translational transition temperatures, i.e. 
all the elements of the theoretical nonergodicity matrix F(q, m) simultaneously jump from 
zero to a nonzero value at the same temperature, or in other words: all degrees of freedom 
freeze at a single temperature. 

In the dipole approximation, i.e. relaxing the diagonality approximation, the equations 
for F(q, m) have been solved in a similar way. In this approximation the value of the critical 
F(q, m) was evaluated stopping the iterative calculation when the average difference between 
two consecutive iterations was of order 10~ 8 . In the dipole approximation, the theoretical 
critical temperature has been found to be about T = 208 K, which is not so different 

from T^f MCT in the diagonal approximation and from the numerical one. As discussed in 
Ref. |L7]|, it is reasonable to suppose that the transition is driven by the diagonal structure 



factors, especially the ones with I = or I = 2, which are the most sensitive to variations 
of temperature. It is thus not surprising to observe that the critical temperature is almost 
insensitive to the introduction of the off-diagonal terms, which display a weaker dependence 
on T. 

The theoretical predictions for F(q,m), in the different approximation schemes, are 
shown in Figs. [1], @, |3| and [| in comparison to the corresponding quantities as evaluated 
from the MD simulation fitting the time evolution of the correlators in the early a-region to 
the von Schweidler law IT7 



S ln , Vn ,(q,m,t) - F, n , JV (g,m) = -H$ Vn ,{q,m)i h + H$ >Vn ,{q,m)t lb + O(P) (28) 

In Fig. p] the theoretical predictions for the COM nonergodicity parameter in the two ap- 
proximations examined in this article are compared to the MD data. Both approximations 
are excellent in the low q range (up to around 30 nm~ l ). The MMCT approximations allow 
the calculation, beside the COM glass form factor, also of the angular nonergodicity parame- 
ters. Of course, in the diagonal dipole approximation only the diagonal elements i^cyo of the 
nonergodicity matrix F can be evaluated. The results are shown in Figs. |2| and |3] together 
with the predictions for the same quantities as evaluated relaxing the diagonality approxi- 
mation. The agreement with the numerical data (symbols in the figures) is satisfactory for 
the I = 1 nonergodicity parameters (both for m = and m — 1), while the predictions 
corresponding to / = 2 are less satisfactory. It is reasonable to expect that the worse results 
obtained for I = 2 are due to the fact that this is exactly the cut-off value, and it is thus 
more sensitive to "boundary" effects. The comparison between theoretical predictions and 
MD data for the off-diagonal terms of F, a comparison which is possible to perform only in 
the dipole approximation scheme, is shown in Fig. [|. Again, the agreement between MD 
data and theoretical results become worse on increasing I towards l co . By comparing the 
results for the two different approximations, we note that both, for the COM (Fig. [I]) and 
angular(Figs. [| and ||) nonergodicity parameters, the coupling to the non diagonal correla- 
tors introduced in the pure dipole approximation contribute very little to the determination 
of the diagonal terms of the glass form factor. This result, as well as the small variation 
of the critical temperature within the two approximations, is due to the small amplitudes 
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of the off-diagonal terms with respect to the diagonal ones and supports the idea that the 
critical behavior of the system is mainly driven by the more intense structure factors. This 
consideration suggests to use only the strongest S[ n! i> n r(q,m) as input of the calculation, 
which may allow to choose higher values for the cut-off l co . 



V. SUMMARY AND CONCLUSIONS 

In the present paper we have performed a quantitative test of MMCT for the SPC/E- 
model for water in the supercooled regime. MMCT is an extension of mode coupling-theory 
for simple liquids to molecular systems. It provides equations of motion for the molecular 
correlators S K)K >(q, t) which form a complete set for any time-dependent two-point correlator. 
Of course, there are infinitely many of them which demands for a truncation of the set of 
MMCT-equations at a cut-off value l co for I and I'. In our case we have chosen l co = 2. 
From a pragmatic point of view this may be justified by the fact that up to today there 
seems to be no experimental method which allows to determine those correlators for / and 
I' larger than two, although this can be done for any numerical simulation. Instead of using 



the molecular correlators one also could use correlators in a site-site representation J20 
MCT for molecular systems in a site-site description has recently been worked out [ITU 
This type of approach has the advantage that for molecules with a few atoms the number of 
correlators is small, e.g. for water there are six correlators in maximum, where one has to 
take into account that both protons can be distinguished for a classical system. A site-site 
description, however, has the disadvantage that it is incomplete, i.e. information has been 
lost. For instance light scattering, dielectric spectroscopy, NMR, ESR, etc. directly yield 
molecular correlators. Whereas the site-site correlators can be represented by the molecular 



ones, the reverse is not possible ||26|| . Hence, from a site-site description no exact information 
can be deduced for the experimental techniques mentioned above. Nevertheless, it would be 
worthwhile to solve the molecular MCT-equations in the site-site representation, e.g. for 
water, in order to compare the critical temperature T c , the critical nonergodicity parameters, 
etc. with the corresponding quantities from MMCT. 

Despite the cut-off for / and /' the MMCT-equations are still rather involved. There- 
fore we decided to perform in a first step a further approximation which is that n and n! 
is put to zero for the static correlators and the nonergodicity parameters. Because of this 
approximation we neglect rotations of both protons around the molecular symmetry axis. 
Within these approximations we have calculated the (unnormalized) critical nonergodicity 
parameters Fi j/ Q (q,m) in the g-frame by solving the corresponding MMCT-equations. As 
we have found that the diagonal, static correlators Sio^o(q,m), compared with the nondiag- 
onal ones, are most dominant with respect to their intensity and temperature-dependence 



17], we have additionally solved the MMCT-equations by assuming all static correlators 
and the nonergodicity parameters to be diagonal in I and The reader should note that 
the diagonality of the nonergodicity parameters is an additional approximation. This lat- 
ter approximation has been motivated by a similar investigation for a supercooled liquid of 



diatomic molecules 22 



The solution of the MMCT-equations yield the critical temperature j v ^ IMCT ^ 208A". 
In the diagonalization approximation we also obtained a single transition temperature 
t mmct 206fC which does not differ much from T^ IMCT and T„ MD = 200 ± 3K. We 
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consider this very good agreement as fortuitous since usually the mode coupling theory 



strongly overestimates the freezing [p4|, P5|. That the diagonalization approximation almost 
does not affect the transition temperature is quite different to what has been found from a 
similar study for diatomic molecules with Lennard- Jones interactions. There the diagonal- 
ization approximation causes the separate freezing of the COM-correlators (I — I' — 0) and 



the "orient at ional" correlators with I = V ^ [p2 |. This quite different behavior is probably 



related to the much stronger translational-orientational coupling in case of water. 

Comparison of the MMCT-results for Fi ji (q,m) with the corresponding MD-results 
leads to the following main conclusions: 

1- -^o,z'o(QS m ) obtained from MMCT, without and with diagonalization approximation, 
differ only slightly from each other which confirms the dominance of the diagonal 
correlators. 

2. The qualitative g-dependence of Fio } i>o(q,m) from the MD-simulation is well repro- 
duced by the corresponding MMCT-result. This is even true for some fine structure 
like the shoulder at q = 20nm~ 1 of F 10iW (q,m = 0). 

3. A good quantitative agreement between the MD- and MMCT-results is found for the 
full g-range for F w ^ w (q, m—1) and Fx 0j00 (q, m = 0). For the other cases, except those 
with I = I' = 2, a reasonable agreement is found for q < 30nm~ 1 . For -^20,20(55^) the 
deviations are rather large, particularly for q < 40nm _1 . 



This behavior for -F]o,z'o(g, m ) is in full accordance with that for diatomic molecules |22| , |25 
The larger discrepancy for the case I — I' — 2 is probably due to the cut-off at l co = 2. For a 
single dumbbell in an isotropic liquid of hard spheres it has been shown that MMCT yields 
accurate results for e.g. Fi 0t i (q, m) if one chooses l co = I + 2 |27] . 



Since this first quantitative comparison of MMCT and the MD-results for water is en- 
couraging we are planning to extend our MMCT-study for water to include n and n'. This 
will offer the possibility to answer the interesting question whether the 180°-jumps of the 
protons, which leave the molecule invariant, freeze at the same critical temperature where all 
the other degrees of freedom freeze, or if they will freeze at a lower temperature. Besides the 
nonergodicity parameters it would also be interesting to calculate from MMCT the critical 
amplitudes Hi n ^ n /(q,m) and the exponent parameter A from which the critical exponent a, 
the von Schweidler exponent b as well as 7 which characterizes the power law divergence 
of the a-relaxation timescale, can be obtained. Finally a solution of the time-dependent 
MMCT-equations would be desirable. This extensions will be hard to achieve without fur- 
ther approximation schemes. Our results give a strong indication that this could be possible. 
As the comparison of the different approximations presented in this paper demonstrates, the 
solutions of the MMCT are not strongly affected, if small and temperature insensitive com- 
ponents of the static structure factors Si n) ii n i (q, m) are neglected. Therefore it should be 
possible to restrict the MMCT to the most relevant components, where the question of rel- 
evance is decided on the basis of the static correlation functions. Depending on the system, 
this procedure can lead to a dramatic reduction of memory functions to be calculated. This 
could enable us to capture more qualitative features of the nonergodicity parameters and 
dynamics, by including the relevant components with / > 2, thereby neglecting irrelevant 
components with I < 2. 
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To summarize, we can say that this first quantitative test of MMCT for water has 
demonstrated reasonable agreement of both, the critical temperature and critical nonergod- 
icity parameters obtained from MMCT and a MD-simulation, although the agreement of 
the transition temperatures should not be overestimated. 
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APPENDIX A: SIMPLIFICATION OF THE MMCT-EQUATIONS FOR WATER 

MOLECULES 

The C^-symmetry of the water molecule leads to the consequence that the molecular 
correlation function S K>K '(q, t) for n and n 1 odd are given by the self-correlations and vanish 
for n odd, n! even or n even,n' odd, i.e. the matrix S(<f, t) has the block structure 





n' even 


n' odd 






S(g,t) 





n even 


S(q,t) = 















n odd 



(Al) 



Thus the direct correlation functions c K ^(q) are nonzero for n and n' even, only. From 
equation (|12|) we can conclude that the functions (q, q\, <f 2 ) vanish if rii and n 2 are 

odd. As a consequence the memory functions can not contain terms of the kind v * v * §( s > * 
S( s \ Using that the factors u (see eqs. flT2]), (|H|)) contain the Clebsch-Gordan coefficients 
C(ZiZ 3 Z; riin 3 n) or C{l\lzl; rii + fi n 3 n + /i) or the corresponding ones with (1 <->• 2), where n 3 
has to be even since it occurs also as index of the direct correlation function (see eq. (|T2]) ) 
one concludes that terms of the kind 

v * v * S * S (A2) 

are only contained in memory functions with n and n' even. Further one finds, that memory 
functions with n and n' odd contain only terms of the kind 

v * v * S * S {s) (A3) 

while memory functions with different n even,n' odd or n odd, n' even vanish. We can 
summarize those findings in the schematic representation 





n! even 


n' odd 






vvSS 





n even 


m(q, t) = 













vvSS^ 


n odd 



Since this block structure is also preserved under matrix inversion the whole set of equations 
is split into two parts. The first part which consists of the block with even n and n' forms a 
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closed set of equation since also in the calculation of the memory kernels (|A2|) only matrix 
elements with even n and n' occur. The second block depends on the solution of the first 
set because of the structure of the memory functions ( |A3| ) and can be shown to be identical 
to the tagged particle equations for the self part. 



APPENDIX B: Q-FRAME REPRESENTATION 

In contrast to the case of simple liquids the density correlator S(g, t) as defined in sec- 
tion U depends on modulus and orientation of the vector q. Therefore a direct numerical 
implementation of the equations given above is not suitable but further reformulations are 
necessary. The dependence on the direction of the "outer" g-vector is easily removed by 
choosing the 2-axis of the laboratory frame of reference in direction of q. The choice of the 
g-frame further offers the advantage that the matrix S(q,t) becomes diagonal with respect 
to m and m', i.e. it holds 

Slmn,l'rn'n'{q&z,t) = S^ lln ,(q,t)5 m>m i = Si n y n > (q, 171, t)5 m%m t . (Bl) 

To get a reformulation of the equations of motion in terms of the g-frame quantities we 
still have to care about the "inner" g^vectors q\ and qi appearing in the summation of the 



MCT memory functions (cf. eq. ([TT]) ■ This can be done by using the transformation law 
of the tensors S(q,t) under rotations which gives a relation between Si n j' n '{q,m,t) and the 
molecular correlation function for general direction of q: 

Slmn,l>m'n'(q,t) = ^ ^inm" (^)^m'm» {^q)Sl n ,l'n' (<?, m" , t) , (B2) 
m" 

where Q q denotes the polar angles of the vector q with respect to the laboratory frame. 
Thus we get the following set of MMCT-equations in the q-frame representation: 

S(g, m, = — zS^ 1 (q, m) + S _1 (g, m)~K(q, m, z)S^ 1 (q, m) (B3) 
Kin,i>n>(q, m,z)=J2Yl qZ(q^z)qv'^'*(qe z )kZ^iCn'(q^ z ) ( B4 ) 



zJ 



-1 



m(q,z) 



(B5) 



with J as given in section |II A| . The memory functions are still nondiagonal with respect to 
m and m! . But one can show that they are different from zero only in the following cases 



a —T 


a' =T 


m + fi = m' + fi' 


a —T 


a' =R 


m + n — vnl 


a —R 


a' =T 


m = m! + // 


a —R 


a' =R 


171 = 171' 



(B6) 



Besides the g-frame representation we have used that the restricted summation over g! and 
g*2, which becomes an integration in the thermodynamic limit V, N — > oo with po — const, 
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can be reduced to a double integral. The general expression for the memory functions in 
the g-frame is given by 



m lmn,Vm'n'Vli 



PO 



In 2 ) 3 Jo 



g+qi 



Q1Q2 



dq 2 — 

k-9i| q 



mim 2 hlfol'z mn'^n'z 



x St^i^q^m^t) Si 2n2) v 2n > 2 (q2,m 2 ,t) 
with the vertex functions 

, B 1 ni,j 3 n 3 (Mi?2;mm 1 m 2 ) = £ [%^ 3n3 ,z 2 „ 2 {qqiq^ mm x m 2 ) c hn3 ,, ini (gi, mi)+ 



+ (- 1 ) m WM 3 n3Ani(?92gi;™2mi) C hm j 2n2 (q 2 , m 2 ) 



The coefficients it are given by 



i 1 ni,i 3 n 3 (95 , i92;mmim 2 ) = l' 1+ ' 2 '(-1) 



(2/! + l)(2/ 2 + l) 



2/ + 1 



<C(«0 x 



x E (^i)^- m - m2 (-OzWtU; m" m - m" m) 



C(hl 2 l; n x n 2 n) 



(B7) 



:bs) 



a = T 

C(lxl 2 l;rii + \in 2 n + fi) a = R 



• (B9) 



The quantities d mn are the reduced Wigner matrices (19| and the angles 9i are determined 
by 



coslfi 



q 2 + q 2 - g\ 
2qq x 



sin(6/i 



1 - cosffl/ 



(BIO) 



and the corresponding relations with (1 2). Further we have 



?te ( 



^ sm(^ 
g i cos(6' i ) 



a = T, n = ±1 
a = T, /i = 



(Bll) 



+ 1) -n(n + fi) a = i?,/i±l 
n a = R, // = 



It should be immediately obvious that a numerical solution of the equations of mo- 
tion given here poses a formidable task. Due to the large number of summations and the 
occurrence of special functions the evaluation of the memory functions will be the main 
computational problem. 

The first step in an analysis of the equations for S(q, m, t) is the localization of the critical 
temperature T c at which a bifurcation of the long-time behavior of the solutions takes place. 
Therefore one studies the nonergodicity parameters 



F(g, m) = hm S(g, m, t) — — lim zS(q, m, z) 



(B12) 



which obey the following equations 
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F(g, m) = [S _1 (g, m) + S _1 (g, m)K(g, m)S _1 (g, m)j 1 



Q/i. 

m lmn,l 



x u Bini,i 2 na (^1^5 mm^) vfy'j,^^ (qqiq 2 ; m'm 1 m 2 ) 
x ^1,^(51,7711,*) F l2n2 j 2n , 2 {q 2 ,m 2 ,t) 

that have to be solved self-consistently. 
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FIGURES 



FIG. 1. Critical center of mass nonergodicity parameter -Foo,oo(<Zj m = 0) as calculated from 
the MD simulation (symbols) compared with the theoretical predictions obtained in two different 
approximation schemes: MMCT in the diagonal-dipole approximation (dashed line) and MMCT 
in the pure dipole approximation (dotted line). 

FIG. 2. Diagonal critical nonergodicity parameters with I = I' = 1 (-Fio,io(#> m)) as calculated 
from the MD simulation (symbols) compared with the theoretical predictions of MMCT in the 
diagonal-dipole approximation (dashed line) and MMCT in the pure dipole approximation (dotted 
line) 

FIG. 3. Diagonal critical nonergodicity parameters with I = I' = 2 (i*20,2o(<7) m )) as calculated 
from the MD simulation (symbols) compared with the theoretical predictions of MMCT in the 
diagonal-dipole approximation (dashed line) and MMCT in the pure dipole approximation (dotted 
line) 

FIG. 4. Off-diagonal critical nonergodicity parameters (I / I') as calculated from the MD 
simulation (symbols) compared with the theoretical predictions of MMCT in the pure dipole ap- 
proximation (dotted line) 
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FIG. 1. L. Fabbian et al 



19 




40 80 

q (nm ) 



FIG. 2. L. Fabbian et al 
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FIG. 3. L. Fabbian et al. 
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FIG. 4. L. Fabbian et al 
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